Spatial Bayesian Distributed Lag Non-linear Modelling of Temperature-Related Mortality in Birmingham: Projections under RCP Scenarios
1 Introduction
The aim of this study was to analyse the association between air temperature and all-cause mortality particularly for Birmingham and to demonstrate the application of spatial Bayesian distributed lag non-linear models (SB-DLNMs) for estimating temperature-mortality relationships at fine geographical level with limited data. Using Birmingham as a case study, this framework can be replicated by other local authorities or cities, both in UK and internationally, to guide their local public health policies.
2 Methods
2.1 Data
Daily all-cause mortality data in Birmingham were extracted from the death register table that was stored in the Birmingham and Solihull Aristotle data warehouse. The dataset covered the period between 01/01/2014 and 31/12/2024 and and included the electoral ward of residence for each deceased individual.To restrict the analysis to Birmingham residents, only those associated with a valid Birmingham ward code were included in the analysis. To avoid bias in the assessment of temperature-related mortality, individuals whose underlying cause of death was COVID-19 (ICD-10 codes U07.1, U07.2) were removed to account for the inflated deaths during the pandemic period. Therefore, a total of 94,684 deaths during the study period were used for analysis..
During the same study period, the daily maximum and minimum values of air temperature on 1 km grid across UK were obtained from the HadUK-Grid database, produced by the UK Meteorological Office. The UK-wide grid data were then spatially subset to Birmingham by cropping and masking the raster layer using the dissolved polygon to exclude all cells outside Birmingham. Daily ward-level minimum and maximum air temperature were calculated using the area-weighted average of the grid cells intersecting each ward polygon. Finally, the daily ward-level mean values of air temperature were derived from the average of the daily ward-level minimum and maximum air temperature.
To capture the future effects of excess mortality under different climate change scenarios, the bias-corrected future series of daily mean air temperature on 1 km grid across UK were obtained from NERC EDS Centre for Environmental Data Analysis, available on CEDA archive [1]. This dataset contains 4 ensemble members (01, 04, 06 and 15) representing different realisations of future climate and under each member, there are four representative concentration pathways (RCPs) namely RCP 2.6, RCP 4.5, RCP 6.0 and RCP 8.5 [2]. These four scenarios correspond to the pathways that describe the future greenhouse gas concentration, where RCP 2.6 represents the best-case scenario with a lower greenhouse gas emission and lower increase in future temperature by 1.6ºC. RCP 4.5 and RCP 6.0 represent intermediate level of greenhouse gas emission with mild increase in temperature by 2.5ºC and 2.9ºC respectively. RCP 8.5 represents the worst-case scenario with a higher greenhouse gas emission and higher increase in future temperature by 4.3ºC. The same procedure was used to spatially subset the projected daily mean air temperature at Birmingham ward level.
2.2 Modelling approach
This analysis adopted the extension of the Bayesian distributed lag non-linear models (DLNMs) to the spatial Bayesian DLNMs (SB-DLNMs) as demonstrated in recent methodological study [3]. For each day \(d\) and ward \(w\) in Birmingham, let \(Y_{d,w}\) denote the observed number of deaths for all causes. The association between ambient air temperature and all-cause mortality was modelled using a Poisson likelihood with a log-link function:
\[ Y_{d,w} \sim \text{Poisson}(\lambda_{d,w}) \]
\[ \log(\lambda_{d,w})=\alpha_{s(d,w)}+\sum_{c=1}^{n_c} \beta_{w,c}\, CB_{d,w,c} \]
where \(\alpha_{s(d,w)}\) denotes the stratum-specific intercepts under the use of case-crossover design, which is a methodology framework extensively used in environmental epidemiology [4]. In order to control for long-term trends, seasonality, and ward-specific time-invariant confounding, the data were stratified by ward, year, month and day of week, such that daily temperature and all-cause mortality in a given ward were compared exclusively with observations from the same ward on the same day of the week, within the same month and year. It has to be noted that strata with zero observed deaths were excluded from the model because they do not contribute to the likelihood and provide no information for parameter estimation; \(CB_{d,w,c}\) is the cross-basis matrix that captures the non-linear relationship between temperature and mortality across a range of lags. Two separate spline basis functions were specified for the temperature and lag dimensions. With regards to the exposure-response space, cubic B-splines were selected with internal knots placing at the 10th, 75th and the 90th percentiles of the temperature range. For the lag-response space, natural cubic splines were used with three internal knots placed at equally spaced values on the log scale, including an intercept. A lag window of 0-21 days was chosen to capture the delayed effect of cold and the short-term mortality displacement from the heat effect. The specification for both basis functions followed previously established practice in the literature [5].
\[ \beta_{w,c} = \beta_c + \beta^{*}_{w,c} \]
\[ \beta^{*}_{w,c}=\frac{1}{\sqrt{\tau_c}}\left(\sqrt{1-\phi_c}\, v_{w,c}+\sqrt{\phi_c}\, u_{w,c}\right) \]
Modelling for a single local authority at fine geographical level often results in extremely high variance and unreliable estimates. To mitigate this, it is assumed that neighbouring wards share similar characteristics that contribute to their underlying mortality risk. The spatial dependency was modelled within a Bayesian hierarchical framework using the Besag–York–Mollié 2 model (BYM2), by incorporating the spatial structure directly into the ward-specific cross-basis coefficients \(\beta_{w,c}\) of the cross-basis matrix \(CB\). To define spatial neighbourhood, queen contiguity was specified such that wards were considered as adjacent if they shared either a long border or just a single corner. Specifically, each coefficient \(\beta _{w,c}\) is modelled as the sum of the global fixed effect (\(\beta_{c}\)) for the entire region, and a ward-specific random effect \(\beta^{*}_{w,c}\) which accounts for the spatial dependence between neighbouring wards. The spatial neighbourhood structure under BYM2 model can be separated into scaled structured (\(u_{w,c}\) ) and unstructured (\(v_{2,c}\)) spatial components. The proportional contribution of these components is controlled by a mixing parameter \(\phi_{c}\), and the overall variability of the spatial effect is controlled by a precision parameter \(\tau_{c}\) [6].
2.3 Prior specification
Since this analysis was conducted within a Bayesian framework, prior distributions of all model parameters were specified as follows:
\[ \alpha_{s(d,w)} \sim \mathcal{N}(0,\; 10^{4}) \] \[ \beta_c \sim \mathcal{N}(0,\; 10^6) \]
\[ \Pr\!\left(\frac{1}{\sqrt{\tau_c}} > 1\right) = 0.01 \]
\[ \Pr\!\left(\phi_c < 0.5\right) = \tfrac{2}{3},\qquad0 \le \phi_c \le 1 \]
It is assumed that the strata specific intercept (\(\alpha_{s(d,w)}\)) and the global fixed effect (\(\beta_{c}\)) follow the weakly informative Gaussian priors. The prior specification for the scaled structured and unstructured spatial components of the BYM2 model uses the Penalised Complexity (PC) priors which were introduced as a interpretable and meaningful specification [7]. For the precision parameter \(\tau_{c}\) , the PC prior was set such that \(\Pr\!\left(\frac{1}{\sqrt{\tau_c}} > 1\right) = 0.01\), implying that substantial there is only 1% prior probability of large residual variation in temperature-mortality risk across wards. For the mixing parameter \(\phi_{c}\), the PC prior was set such that \(\Pr\!\left(\phi_c < 0.5\right) = \tfrac{2}{3}\), reflecting a prior specification that there is a 67% probability of less than half of the residual variation in temperature-mortality risk across wards is spatially structured. This means that unstructured, ward-specific random noise accounts for a greater proportion of the variability. These conservative priors were adopted because BYM2 model naturally shrinks towards \(\phi_{c}=0\) unless the data provide strong evidence of genuine spatial dependency [8].
2.4 Statistical analysis
The associations between temperature and mortality across all wards were estimated with a single-stage approach with the SB-DLNM. The model was fitted using the R-INLA (Integrated Nested Laplace Approximation) package which provides computational efficiency in execution times than using Markov Chained Monte Carlo method from other R packages. Given the high dimensionality of the model induced by the cross-basis matrix and ward-specific spatial effects, inference was conducted using empirical Bayes strategy instead of the full Laplace approximation.
The output of the SB-DLNMs were obtained by drawing 1,000 independent samples from the posterior distributions of the model parameters for each ward. These posterior samples were then used to derive the ward-specific temperature-mortality associations by aggregating the coefficients of the global fixed effect for the entire region and the ward-specific cross-basis coefficients. In order to calculate ward-specific relative risk (RR) curves for each ward, grids of ward-specific temperatures were constructed using the observed minimum and maximum temperatures within each ward. These grids were then used to generate prediction cross-basis matrices for each ward. These matrices were then combined with the corresponding prediction cross-basis matrices for each ward to reconstruct the cumulative temperature–mortality association across the observed temperature range. Initially, the reference point of the temperature was selected at the 90th percentile which is around 17 to 18 °C that matches the nationwide reference under UK context.
To obtain the minimum mortality temperature (MMT) and the MMT distribution, the 1,000 RR curves calculated from the previous step were re-centred by finding their minimum RR for each ward. The search for the minimum RR on each curve was constrained within the 1st and the 99th percentiles of the ward-specific temperature grid. The constraint was justified because extreme temperatures at extreme tails of the distribution are rare, excluding them can provide more stable estimates, which is in line with previous studies [9]. Once the minimum RR values were identified, the corresponding temperatures from the temperature grids were extracted and defined as the MMT. The posterior distribution of the MMT for each ward was the obtained by identifying the MMT across all posterior samples.
Once ward-specific MMTs had been identified for each posterior RR curve, excess mortality attributable to heat and cold were then calculated using the backward attributable fraction (b-AF) framework [10]. The posterior distributions of the cumulative temperature-mortality association across lags that were centred at ward-specific MMTs were reconstructed by combing the posterior samples of the cross basis coefficients. The cumulative log-relative risks were calculated by summing lag-specific contributions over the 21 lags and transformed into b-AF for each posterior draw for each ward. These b-AF(s) were the multiplied by the daily observed deaths to derive the attributable numbers (ANs) for each posterior draw for each ward. Daily temperatures that were above the MMT(s) were counted as cold-related and above as heat-related.The cold- and heat-related excess mortality were calculated by aggregating the ANs separately for cold and heat and by year. Finally, the annual cold- and heat-related excess mortality were obtained by averaging over the study period, and standardised to per 100,000 in population using the 2022 ward level population estimates.
For projections, 1,000 posterior distribution of the coefficients of the ward-specific temperature-mortality associations from the fitted model were extracted along with the ward-specific baseline daily mortality for each ward obtained by refitting the model to the observed series and averaging the resulting expected deaths. Since there ward-level population projection were not available, this analysis did not model future demographic change or adaptation, implying that the the results should be interpreted as the potential impacts of climate-driven change under a constant population and vulnerability assumption rather than predictions. The baseline impacts were calculated over 2014-2024 (using the same calender window as the study period) and extracted from each climate scenario. Future impacts for each scenario were then estimated for the 2030s, 2040s, 2050s, 2060s and 2070s decades. The procedure to obtain the future ANs and standardised excess mortality was the same as described in the previous paragraph. Given that each climate scenario projection consists of 4 ensemble members, the results were pooled posterior draw across all ensemble members.
Uncertainty was quantified using posterior median and 95% credible interval (CrI) obtained from the posterior distributions. All data processing was performed in R version 4.4.2. All statistical analyses were performed in R version 4.4.2. All data visualisations were performed in R version 4.4.2. The code is available on ….
3 Results
Figure 1 shows the overall ward-specific cumulative relative risk curves for an adjacent set of neighbouring wards in Birmingham. It is observed that these RR curves show the usual J-shape with varying MMT but generally occur at a relatively warmer temperatures at around the 90th percentile temperatures. Elevated mortality risks in both cold and heat temperatures are observed across wards.
Figure 2 shows the spatial distribution of the median RR at 1st and 99th percentile of the ward specific temperatures referenced to their own MMT. Billesley exhibits the highest RR at 1st percentile temperature at 2.58 (95% CrI: 1.47-4.61). In contrast, the RR values at 99th percentile are also elevated where ward Glebe Farm & Tile Cross shows the highest RR at 1.83 (95% CrI: 0.99–3.24). The cold-and heat- extreme map shows consistent elevated RR across wards, with similar RR values between adjacent neighbours, implying spatial clustering and vulnerability to lower and higher temperatures is not uniform across Birmingham.
Figure 3 shows the posterior median annual cold and heat-related excess mortality per 100,000 population. It shows similar pattern where areas with higher excess mortality also have higher RR being reported in figure 2. Figure 4 shows the annual excess mortality per 100,000 for temperatures below 2.5th and above 97.5th percentiles, it demonstrates completely different pattern especially for hotter temperatures where elevated excess mortality is observed in Bilesley. For colder temperatures, the pattern is slightly different where Castle Vale dominates the excess mortality.
Figure 5 and 6 show the cold- and heat-related excess mortality at ward-level projected under 4 different RCPs . The colour is getting lighter for cold-related excess mortality and conversely getting darker for heat-related excess mortality. The projected change under RCP 8.5 shows the strongest shifts with clearer decline in cold-related excess mortality and heat-related excess mortality by the 2060s and 2070s. Under lower emission pathways (RCP2.6-6.0), changes are less obvious and smaller. Since the projection used the fitted coefficients from the model, the spatial pattern remains similar to previous figures. Overall, these figures illustrate a transition from cold-dominated burden toward increasing heat-related risk later in the century.
4 Discussion
This study adopted the novel SB-DLNM for estimating the short-term health impact of non-optimal temperature for Birmingham at fine geographical level and obtained more stable estimates with the use of a single stage model that take into accounts of spatial dependences. Uncertainties associated with model parameters and epidemiological metrics were fully propagated within the Bayesian framework. This study also attempted to project the future excess mortality under different scenarios to provide a comprehensive assessment to inform and shape policies at local authority level.
The results are broadly comparable with previous studies focusing on nation-wide temperature-mortality association. Almost all of the wards in Birmingham were estimated with MMT(s) that are very close to the range of “comfortable” temperatures that suggested and studiied in UK. [5] However, small numbers of wards exhibited lower MMT(s) which can lead to extremely high heat-related excess mortality. This is because any temperature above MMT are calssified as heat-related and therefore wards with lower MMT(s) will classify a larger proportion of days as “heat”, increasing the heat-attributable deaths. Consequently, the plot that showed excess mortality below temperatures at 2.5th and above 97.5th percentile displayed completely different patterns. It has to be noted that under the UK climate, the 2.5th and 97.5th percentile of daily mean temperatures correspond to relatibely common hotter and colder days. rather than extreme conditions.
The projection is in line with previous studies where decreasing cold-related mortality and increasing heat-related mortality are estimated for the future. Under the most extreme RCP 8.5 scenario, the shift is most pronounced where warming is greatest and the divergence between cold- and heat-related excess mortality becomes clearer. As mentioned previously, the projected results did not consider adaptation of population change and therefore results should not be interpreted as predictions. If adaptation were considered, the projected result might draw different conclusion. Recent European studies that incoperate population change suggest that cold-related excess morality can increase before the first half of the century despite warming, it is because aging population increases the pool of vulnerable individuals and increases baseline mortality in older age groups. Given that the subnational population projection for Birmingham also indicate aging population, cold-related excess mortality in Birmingham might plausibly increase in the near future even if temperature is getting warmer.
This study has some limitations and the audiences must interpret with caution. Firstly, the estimated CrI(s) were extremely wide for both the observed-period estimates and the projected results, reflecting uncertainty in temperature-mortality relationship with limited data. Fortunately with the use of SB-DLNM, posterior medians are comparatively stable. As a result, the estimates should be interpreted as plausible ranges and broad patterns which are useful for identifying spatial heterogeneity rather than precise estimates for individual wards. Secondly, this study did not use directly age-standardisation but instead the unadjusted crude rate which might underestimate or overestimate the excess mortailty rate. The choice was driven by the following two constraints. Ward level population projection are not available and the further separate estimation of temperature-mortality relationship for each age group would substantially increase data sparsity, leading to more unstable estimates. Thirdly, projection with no adaptation held the baseline mortality constant, this isolate purely the impact of climate change but less meaningful to guide immediate policy making given that the population in Birmingham is likely to age. Finally, this study only used 4 out of 12 ensemble members of the UK climate projection which might not capture parts of the internal variability and might not represent the full range of plausible future climates.
5 Conclusion
In conclusion, the study provides a practical example of utilitsing SB-DLNM with limited data to obtain at least a broad picture of the temperature-related mortality within a local authority with fine geographical level. By producing ward-level estimates and projections with quantified uncertainty, the approach offers a “broad” idea of where the heat and cold burdens most likely to be concentrated, supporting tailored policy making. More importantly, other local authority in the UK can also implement the same method to guide their local policy without relying on large scale national study which updates rarely. Finally, it should also be noted that local authorities that have access to deaths data further back in years will be able to obtain the temperature-mortality associations with greater precision, leading to more robust estimates and greater certainty.
References
[1]
Robinson EL, Huntingford C, Semeena VS, Bullock JM. CHESS-SCAPE: Future projections of meteorological variables at 1 km resolution for the United Kingdom 1980–2080 derived from UK Climate Projections 2018 2022. https://doi.org/10.5285/8194b416cbee482b89e0dfbe17c5786c.
[2]
IPCC. Summary for policymakers. In: Pörtner H-O, Roberts DC, Masson-Delmotte V, Zhai P, Tignor M, Poloczanska E, et al., editors. IPCC special report on the ocean and cryosphere in a changing climate, Geneva, Switzerland: IPCC; 2019.
[3]
Quijal-Zamorano M, Martinez-Beneito MA, Ballester J, Marí-Dell’Olmo M. Spatial bayesian distributed lag non-linear models (SB-DLNM) for small-area exposure-lag-response epidemiological modelling. International Journal of Epidemiology 2024;53:dyae061. https://doi.org/10.1093/ije/dyae061.
[4]
Tobias A, Kim Y, Madaniyazi L. Time-stratified case-crossover studies for aggregated data in environmental epidemiology: A tutorial. International Journal of Epidemiology 2024;53:dyae020. https://doi.org/10.1093/ije/dyae020.
[5]
Gasparrini A, Masselot P, Scortichini M, Schneider R, Mistry MN, Sera F, et al. Small-area assessment of temperature-related mortality risks in england and wales: A case time series analysis. The Lancet Planetary Health 2022;6:e557–64. https://doi.org/10.1016/S2542-5196(22)00138-3.
[6]
Moraga P. Spatial statistics for data science: Theory and practice with r. Boca Raton, FL: Chapman & Hall/CRC; 2023.
[7]
Simpson DP, Rue H, Martins TG, Riebler A, Sørbye SH. Penalising model component complexity: A principled, practical approach to constructing priors 2015.
[8]
Riebler A, Sørbye SH, Simpson D, Rue H. An intuitive bayesian spatial model for disease mapping that accounts for scaling 2016.
[9]
Tobı́as A, Armstrong B, Gasparrini A. Investigating uncertainty in the minimum mortality temperature: Methods and application to 52 spanish cities. Epidemiology 2017;28:72–6. https://doi.org/10.1097/EDE.0000000000000567.
[10]
Gasparrini A, Leone M. Attributable risk from distributed lag models. BMC Medical Research Methodology 2014;14:55. https://doi.org/10.1186/1471-2288-14-55.